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ABSTRACT 

in 

\ A new code and methodology are introduced for solving the general relativistic magneto- 

ns) ' hydro dynamic (GRMHD) equations in fixed background spacetimes using time-explicit, finite- 

pi ■ volume discretization. The code has options for solving the GRMHD equations using traditional 

artificial-viscosity (AV) or non-oscillatory central difference (NOCD) methods, or a new extended 
AV (eAV) scheme using artificial-viscosity together with a dual energy-flux-conserving formula- 
tion. The dual energy approach allows for accurate modeling of highly relativistic flows at boost 
factors well beyond what has been achieved to date by standard artificial viscosity methods. It 
J> ■ provides the benefit of Godunov methods in capturing high Lorentz boosted flows but without 

complicated Riemann solvers, and the advantages of traditional artificial viscosity methods in 
their speed and flexibility. Additionally, the GRMHD equations are solved on an unstructured 
■ grid that supports local adaptive mesh refinement using a fully threaded oct-tree (in three dimen- 

sions) network to traverse the grid hierarchy across levels and immediate neighbors. A number 
of tests are presented to demonstrate robustness of the numerical algorithms and adaptive mesh 
framework over a wide spectrum of problems, boosts, and astrophysical applications, including 
relativistic shock tubes, shock collisions, magnetosonic shocks, Alfven wave propagation, blast 
waves, magnetized Bondi flow, and the magneto-rotational instability in Kerr black hole space- 
times. 
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!_i ' 1. Introduction 

. P.. 

Magnetohydrodynamic (MHD) driven processes in high-energy astrophysical environments can be ex- 
tremely difficult to model numerically, especially where flows are strongly trans-sonic, characteristic speeds 
approach the speed of light, or nonlinear effects couple over short temporal and large spatial dynamical 
ranges. For example, studies of phase transitions in the early universe (Fragile & Anninos 2003), and the 
formation of primordial magnetic fields and black holes demand resolution across many decades of scale to 
model phase interactions, wall dynamics, and microphysical coupling effects. Similarly, simulations of black 
hole accretion flows (e.g. Fragile & Anninos 2005), thought to exist in quasars, active galactic nuclei, X-ray 
binaries, core-collapse supernovae, and gamma-ray bursts, require stably evolving the GRMHD equations 
in a strongly curved background spacetime over many dynamical scales. Addressing such problems requires 
numerical algorithms that are robust enough to simulate ultra-relativistic waves and shocks, and accurate 
enough to treat various physical instabilities that can be seeded at very small amplitudes and frequencies 
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and yet affect global solutions. We also desire algorithms that can capture the propagation of outflows or 
highly compressed features beyond the limits imposed by single mesh simulations. 

A number of codes have been developed in recent years to solve both the special and general relativistic 
MHD equations. Our approach is most similar to that taken by De Villiers & Hawley (2003a), which is 
based on an internal energy formulation with artificial viscosity (AV) for shock capturing and a method 
of characteristics for the magnetic fields. However, other approaches include the simplified total variation 
diminishing (sTVD) method (Koide et al. 1999), and approximate Riemann solver schemes (Komissarov 
1999; Gammie et al. 2003; Duez et al. 2005; Anton et al. 2005) based on fully conservative formulations of 
the MHD equations. Riemann methods are significantly more complex than AV schemes, but are considered 
more robust for ultra-relativistic problems. Although viscosity-based methods are generally accurate, fast, 
simple to implement, and easy to expand to include multi-physics capabilities, they have historically been 
limited in their range of applicability to moderately relativistic flows, with boost factors less than a few. 

In this paper we describe a new massively parallel, multi-dimensional numerical code Cosmos+ + with 
improved (over traditional AV methods) shock capturing capabilities in the high boost regime and adaptive 
mesh refinement (AMR), representing a significant advance over our previous numerical code (Cosmos , 
Anninos & Fragile 2003). Cosmos++ introduces a more complex unstructured mesh system that allows 
for arbitrarily connected hexahedral (quadrilateral in 2D) cells to conform to any boundary or shape. It 
also allows for dentritic-type zoning characteristic of refined grids at level interfaces. Our AMR framework 
differs from the more standard approach (Berger & Oliger 1984; Berger & Colella 1989) by refining indi- 
vidual cells rather than introducing patches of sub-grids composed of multi-dimensional arrays of cells, thus 
providing greater flexibility in modeling complex flows and greater efficiency in positioning computational 
resources. This framework is similar to that described by Khokhlov (1998), but we generalize the method to 
unstructured meshes. Also, as in our predecessor code, several schemes have been implemented to solve the 
hydrodynamical equations, including both artificial viscosity and non-oscillatory central difference (NOCD) 
methods. However, here we use discrete finite-volume methods in place of finite differences due to the un- 
structured nature of the grid. We also introduce in this paper a new dual energy procedure (eAV) that 
allows conventional artificial viscosity methods to be extended and work robustly at arbitrarily high Lorentz 
factors. 

The basic equations, numerical methods, and tests of the code are described in the remaining sections. 
Although Cosmos++ supports many different physics packages in both covariant Newtonian and general 
relativistic systems, including chemical networks, flux-limited radiation diffusion, self-gravity, magnetic fields, 
and radiative cooling, only the GRMHD algorithms in the AMR framework are covered in this paper. The 
Newtonian physics capabilities are currently the same as those of Cosmos , which were presented in Anninos 
ct al. (2003) and Fragile et al. (2005). 

2. Basic Equations 

Two separate formulations of the GRMHD equations are presented in this section. The first form 
incorporates an evolution equation for the internal energy that is appropriate for schemes using artificial 
viscosity methods for capturing shocks (Wilson 1972, 1979; Hawley et al. 1984). The second form is derived 
directly from the primitive form of stress-energy conservation and is thus fully conservative and provides 
the basis for the non-oscillatory central difference schemes. We use the standard notation in which 4(3)- 
dimensional tensor quantities are represented by Greek(Latin) indices. 
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For a perfect fluid, the stress-energy tensor is generated from a linear sum of the hydrodynamic 



and magnetic contributions Tg 

= phu^u v + Pg^ + -L (i^HBH 2 + u f *u u \\B\\ 2 -B»bA (1) 
= (ph + 2P B )uV + (P + P B )g^-^-B»B 1 ' , 

h=1+e+ r + m =1+Te+ m (2) 
p p p 

is the relativistic enthalpy, e is the specific internal energy, P is the fluid pressure, ||Q|| is the bulk scalar 
(artificial) viscosity, Pb — ||i?|| 2 /87r = g^B^B" /Sn is the magnetic pressure, and B^ is the magnetic 
induction in the rest frame of the fluid. Also, is the contravariant 4-velocity, g^ u is the 4-metric, and T 
is the adiabatic index assuming an ideal gas equation of state P = (T — l)pe. 



where 



2.1. Internal Energy Formulation 

This method uses a form of the GRMHD equations similar to De Villicrs & Hawlcy (2003a), derived 
from velocity normalization u^u^ = — 1, baryon conservation V M (pw Ai ) = 0, energy conservation w„V = 
0, momentum conservation (g av + u a u u )\7 = 0, and magnetic induction V AJ (u A1 i? l/ — B^u v ) = 0. 
Nevertheless, it is worth explicitly writing these equations out for clarity, since we employ a unique expansion 
and grouping of the terms in our numerical implementation. In flux-conserving form, the evolution equations 
for mass, internal energy, momentum, and magnetic induction are: 

d t D + fy(DV l ) = 0, (3) 
dtE + diiEV*) = -(P+ k w \\Q\\)dtW- [P5j + Qj) d.iWV 1 ) , (4) 

dtSj+diiSjV*) = ^dtiV^BjB ) + L di{ ^— gBj B l ) 

47T 47T 

__ _ ^B»B") djg^ J-g di((P + P B )5) + Q)) , (5) 

d t B j + diiBW 1 ) = B%V j + v , (6) 

where g is the 4-metric determinant, Sj is the Kronecker delta tensor, W = ^J~^gvP is the relativistic boost 
factor, D = W p is the generalized fluid density, V % — u % /u° is the transport velocity, = W(ph + 2Pb)u^ 
is the covariant momentum density, E — We — Wpe is the generalized internal energy density, Qj is the 
tensor artificial viscosity used for shock capturing, is a switch used to activate a viscosity multiplier for 
the d t W source, and r\ is a coefficient related to the largest characteristic speed in the flow multiplying the 
divergence cleanser function. Notice there are two representations of the magnetic field in these equations: 
B^ is the rest frame magnetic field 4- vector defined in the stress tensor definition (1) and 

B" = W(B^ - BV) (7) 



-4 - 



is the divergence-free {dB l jdx % = 0), spatial [B = 0) representation of the field. The time component of 
the magnetic field B° is recovered from the orthogonality condition B^u^ = 



B u = - 



and we use 
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(8) 
(9) 

(10) 



to compute the magnetic pressure from the divergence-free field. 



We note that the term di(y/—gBjB l ) in equation (5) can be cast into projected and divergence com- 
ponents B l di(yJ—gBj) + y/ —gBjdiB 1 . The advantage of this form is that the method of characteristics 
described by Stone et al. (1992) can be applied easily to the projected component after generalizing the defi- 
nitions of density and momentum accounting for relativistic inertia and boost contributions. This procedure 
is important for achieving stable simulations of sheared Alfven waves, as demonstrated in §4.2.1. 

Two additional sets of equations are needed for the transport velocity V and boost factor W. The 
velocity is derived from the 4-momentum normalization S^S^ = W 2 (ph + 2P B ) 2 u li u^ = —W 2 {ph + 2P B ) 2 , 
from which we compute 



So 



g 0l S t J_ 

~g 6r + ~g™ 



(g 0i Si) 2 - .g 00 (W 2 {ph + 2P B f + gVSiSj) 



1/2 



(11) 



and then V 1 = S t /S°. A convenient formula is derived for the boost factor using u° obtained from the 
4-velocity normalization written as = u V i S ll /(W(ph + 2P B )) = -1. The boost (W = xf-gu°) is 

then evaluated in one of two ways 



W = — 



(-S„S")V2 



(12) 



where V° = 1, S° is computed from the 4-momentum normalization described above, and the other fields 
are known from solutions to the evolution equations. 



2.2. Conservative Energy Formulation 

A second class of numerical methods presented in this paper are based on a conservative hyperbolic 
formulation of the GRMHD equations. It is the same approach used by Anninos & Fragile (2003), except 
here we include magnetic fields. The equations are derived directly from the conservation of stress-energy 
V Al T' il/ = and then decomposed into space and time components 

d t (J—g T 0l/ ) + d^^—g T lv ) = (13) 

with curvature source terms 

S" = -V^9 r^ 7 . (14) 
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The form of the differential equations that we solve, derived after substituting the perfect fluid stress tensor 
(1) into (13), are 



d t s j + d i {s j v i ) + d i F ij = t? , 



(15) 
(16) 



where we have explicitly split off the transport term from the other divergence flux contributions F la , defined 
as 



1 



F la = J~g \{g La - g^V 1 ) (P + P B ) - ^O^ 5 " ~ B°B a V l ) ) . 
The variables D, V 1 , and g are defined as in section §2.1, and 

W 2 



£ = V^gT 00 
s 3 = V=gT°j 



- ,Mph + 2P B ) + g oa (P + Pb) - —V=gB a B° , 

^L(ph + 2P B )V 3 + y/=g g°i(P + P B ) - -^-\^gB° B j 



(17) 



(18) 

(19) 
(20) 



are the new expressions for energy and momentum. We also note that the divergence-free magnetic induction 
equation (6) can be written in the fully conservative form 



d t B j + di{B 3 V l - B l V j ) = j] , 

as required by the central difference schemes described in §3.2. 

ft is convenient to express £ and S l in terms of the internal energy fields E and S a 



(21) 
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The inverse energy relation 
E ( £^/=g D 
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(22) 
(23) 

(24) 



is useful for computing the fluid pressure (P = (T — 1)E/W), and for comparing relative errors in the dual 
energy update procedure described in §3. 



3. Numerical Methods 

Cosmos++ is designed using object-oriented principles to create mathematical abstraction classes for 
vector and tensor (both three and four dimensional) objects on which most functional operations are based. 
Cosmos++ also takes advantage of the operator overload, inheritance, polymorphism and virtual methods 
features of the C++ language in its design of all the basic class structures. This simplifies the user interface 
considerably, reduces the amount of coding, and allows for the code to be easily developed and expanded. 
Figure f illustrates the relative dependencies of the basic classes, ranging from the fundamental math, 
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MPI, zone and mesh classes, to the intermediate structures defining the metric geometry, field operators 
(e.g., gradients, filters, etc.), and boundary conditions, leading finally to the higher level classes where the 
physics, output, refinement, and evolution driver algorithms depend on all the lower level functions. The 
classes outlined with solid double lines in Figure 1 represent base classes where virtual methods functions 
proved advantageous in their polymorphic design. Classes outlined with double dashed lines indicate where 
polymorphism can help simplify and improve the code design as we add new options for background metrics 
and gradient operators. 

All data structures (mesh and fields) are stored as vector or map container classes using the C+ + 
Standard Template Library (STL) which provides convenient memory access functions and other useful 
support features (e.g., sorting, searching) for interacting with the data. Multi-dimensional STL vectors are 
used as storage containers for pointers to zone objects at each refinement level, which include all the field, 
cell and connectivity data as class member attributes. STL maps provide an efficient mechanism to identify 
direct and inverse relations between zone iterators and unique global zone identifications. 

The essential cell geometry used in constructing meshes for Cosmos++ is a quadrilateral (hexahedral) 
shape in two (three) spatial dimensions. In order to provide enough flexibility so individual cells can be 
arranged in an arbitrarily distorted and unstructured fashion, allowing even reduced or enhanced nodal 
connectivities, we store a number of attributes for each cell, including node positions, inward (towards the 
cell center) pointing face area normals, and zone volumes for convenience (see Figure 2). Each zone also 
carries pointers to all of its neighbors sharing a common cell face. Where a neighbor cell has been refined 
to a higher level, the pointer references the neighbor's parent so that neighbors arc always at or below the 
refinement level of the referencing cell and there is at most one neighbor for each cell face. The actual 
operational neighbors are found by selecting the neighbor's appropriate children if they exist. Neighbor 
searches are carried out each time a refinement cycle adjusts the mesh, and performed with local tree scans 
by going up to the parent, across to the neighbors, then down through the child zone lists. This is done 
to allow greater flexibility and simplicity during the refinement cycle in case, for example, a cell and its 
neighbors are tagged for de-refinement simultaneously. When a cell is refined, it is decomposed into two, 
four or eight sub-cells in one, two or three dimensions, and pointers to these newly allocated child zones are 
stored in ordered contiguous fashion for fast and easy reference. Each child zone, in turn, carries a pointer 
to its parent cell. As a result, all cells in the mesh are fully threaded with local binary-, quad-, or oct-tree 
parent-child hierarchies and semi-direct neighbor access. 

Cosmos++ currently supports several refinement criteria, including field value, normalized field gradient, 
normalized field curvature, mass, and Jeans mass. When either of the field, gradient, or curvature options are 
selected, the user can specify any stored variable upon which to apply the criterion. Also, as an option, the 
refinement criterion can be computed for either the evolved fields or their conformal counterparts (by dividing 
out the metric determinant), which is useful to track only flow (not metric) features. Fields in newly created 
refined cells are reconstructed by a linear, conservative, and monotonic interpolation of nearest neighbor and 
parent data. 

Numerical calculations are performed in the physics packages only across the list of leaf zones which 
have no children. The effective mesh is thus unstructured in general, composed of many different sized zones 
containing any number of hanging nodes with reduced nodal connectivity. In order to accomodate these 
kinds of meshes (and even more general ones) we use finite volume, in place of finite difference, methods to 
solve the GRMHD differential equations. Simple derivatives are computed using standard second order (on 
uniform orthogonal meshes) finite volume discretization. Gradients of a generic vector field T t , for example, 
are evaluated by averaging the gradient function over a single cell control volume and converting the volume 
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integral into a surface summation 

~\ f' ~1 C "I faces 

Gij = W = - djTi dV = --j>T i dA j = --Y^ Aj) f , (25) 

where the summation is performed over all cell faces. Here V z is the zone volume, Aj is the area vector normal 
to the cell face (/) pointing inwards toward the cell center, T* is the appropriately averaged (or upwinded) 
field value at the cell faces, and the gradient is returned as a tensor with rows (columns) representing the 
vector (gradient) directions (in usual gradient index notation — T^j). The negative sign in (25) is a 
result of choosing an inward pointing area vector. At zone interfaces separating regions of different refinement 
levels, we use volume weighted averages of fields to define the "downwind" component for interpolation to 
the cell faces. 

Many of the gradient calculations appearing as source terms in a typical time sequence update are 
computed using a variation of the above expression. Some of the other, more specialized, operations used in 
solving the GRMHD equations in the different formalisms are described below. 



3.1. AV Method 

3.1.1. Advection 

Advcction is solved for each evolved field quantity using a time-explicit, first order forward Euler scheme 
together with the scalar divergence form of equation (25). Letting F represent any of the evolved fields (D, 
E, Sj, B J ), the discrete, finite- volume representation of the transport source follows from (25) 

faces 

diFV*) = £ ( f * yi A *h - ( 26 ) 
Vz f 

where V z is the local donor zone volume, (Ai)f is the inward pointing area normal vector of face /, and 
{V 1 ) f is the face-centered velocity defined as a weighted average across neighboring cells. The quantity (F*)/ 
represents first order zone-centered fields estimated at each cell face by a Taylor's expansion using limited 
gradient extrapolants, F* = F z + (diF)^(r 1 —r\), projected from the donor cell center r\ to either the face 
center r % = r l j or the advection control volume r % = r % f — (At/2)(V l )f, over a time-step interval At. 

The zone-centered limited gradient (djF)^ is constrained to force monotonicity in the extrapolated 
fields. In the case of unstructured meshes, this is achieved by identifying three unique control volumes 
which we assign as upstream, downstream, or average. The average control volume is simply the total 
cell volume. As shown in Figure 3, upstream and downstream volumes are constructed from the sub-zonal 
tetrahedral (triangular) in 3D (2D) geometric elements composed of node positions connecting a cell face 
center, two (one) zone nodes in 3D (2D), and the zone center. A sub-zone associated with the cell face / 
is defined as upstream (downstream) if the sign of the vector product (AiV 1 )/ is positive (negative). The 
total upstream and downstream volumes are the sum of the tetrahedral (triangular) sub-zones matching the 
corresponding signature criteria. Gradient operators are then computed with (25) on each of these unique, 
arbitrary polyhedral (polygonal) sub- volumes, as surface integrals with appropriately averaged fields at each 
surface element boundary. We denote these gradients as dfF, dfF, and dfF for the upstream, downstream, 
and average, respectively. To enforce monotonicity, the actual limited gradient (c^F)^ is set to zero if the 
vector product of any combination of the three gradient operators is negative. The gradient is limited further 
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to various degrees of sharpness by defining the normalized scalar 9g — dfF / '(dfF + e), for each direction 
£ and smallness parameter e <C 1, and applying cither of the minmod <\>i = max[0, min(l, 6t)], van Leer 
4>t = (1^1 +^)/(l + \0t\), or supcrbec (f>e = max[min(l, 26 g), min(2, 9g)] limiters (LeVeque 1992), to derive 
the final expression for the zone-centered limited gradient components (dgF)^ = <j)edfF. 

An alternative, though somewhat more restrictive and costly, method that we have developed for 
applying unstructured grid limiters on gradient functions is based on modifying the magnitude of the 
average gradient with some function of the maximum (<9 max = maxfldf 7 !, \df\, \df\]) and minimum 
(9 m ; n = minfl^l, \df\, \df\]) of the three sub-zonal gradient magnitudes. For example, 

(0^ = ah (/3d max , (1 - (3)d min ) J^ F , (27) 

where (3 is a steepness parameter bounded by < /3 < 1, h(...) is any function of the arguments, and a is a co- 
efficient to enforce monotonicity in the extrapolated field Fj. In particular, we set a = min(l, max(0, min(ai, a 2 ))), 
where 

r max (o, f„„-f,) if F/ = Fz + ( a.F)^(r* -rj)/a>ma X (F Z) F Zin ) 



and 



[ 1 otherwise 

1 otherwise 



where F Zjn refers to the field value in the neighbor zone center. The min/max operations in (28) and (29) 
are performed over each of the faces in the donor cell, and a is chosen as the smallest value needed for 
strict monotonicity across all cell faces. Of course these procedures for solving the advection equation and 
evaluating limited gradients reduce in one spatial dimension to the familiar upwind scheme as described, for 
example, in Hawley et al. (1984). 



3.1.2. Artificial Viscosity 

We have implemented five different artificial viscosity options for shock capturing. All of these constructs 
have a common inertia multiplier defined for relativistic MHD as 

I = D + E + W(P+\\Q\\ + 2P B ) . (30) 

The inertia is also normalized by a factor Ip that scales out the local 3-geometry in curvilinear coordinates, 
and introduces a relativistic boost factor to an arbitrary power that is useful for effectively transforming the 
length scale between proper and boosted frames, and providing a flexible scaling with boost. The normalized 
inertia multiplier is thus defined as 

/ ,„ ..... i AO 



In = j^ = (D + E + W(P+\\Q\\ + 2P b )— , (31) 

where 7 is the 3- metric determinant, and n is an arbitrarily specifiable constant. 

The simplest artificial viscosity we consider is the scalar form based on the work of von Neumann & 
Richtmyer (1950): 

Qi ^(l N Ald k V k (k q Ald k V k -k l C s )S^ if 9^<0 
I otherwise 
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where Al is the minimum covariant zone length, C s is the local sound speed, and k q and ki are constant 
coefficients multiplying the quadratic and linear contributions, respectively. This scalar viscosity remains 
one of the most popular representations due to its simplicity and robustness; however, it is also the most 
diffusive of the methods we have implemented since it indescriminantly filters out transverse, longitudinal, 
and sheared isotropically averaged compressive flows. 

A second option for viscosity is similar to (32) but is designed for anisotropic shock capturing along the 
principal grid axes similar to what is achieved in directionally split, fully structured mesh approaches: 

Qj = Iji ^(k q ai o t v l - h c s ) diag&v* - \d x v% a v vy - \d v vy\, d z v z - \d x v*\] , (33) 

where diag[...] refers to a list of the main diagonal elements of a tensor. This is generally less diffusive than 
the isotropic scalar viscosity option since it is automatically disengaged along any grid direction that is not 
undergoing compression. 

Another option is based on the scalar viscosity developed by White (1973), but extended here for 
relativistic problems. Defining the quantity 

Q = Al |^H 1/2 \9 iJ diP ^(V 1 )! 1 / 4 , (34) 

this viscosity can be written as 

Q) =I N Q{k q Q + h C s )5) (35) 

when both conditions, g^diP dj(I^~ ) < and diV 1 < 0, are met. The advantage of (35) is that both pressure 
and velocity gradients are incorporated into its definition, which is a better indication of the presence of 
shocks than the velocity gradient by itself. It can therefore be more effective at suppressing artificial heating 
in regions undergoing adiabatic compression. 

The final two options we consider are genuine tensor viscosities. The first of these is similar to Tschar- 
nuter & Winkler (1979), but does not include the full covariant gradient treatment. However, since the 
velocity gradients are computed in the conformal frame and the proper volume is accounted for in the inertia 
normalization, the version we have implemented 

Q) = I N Al {k q Al d t V l - h C s ) Sym (d^V* - ^<*j-) , (36) 

where Sym(...) denotes a symmetry operation, is a reasonable simplification. This traceless tensor viscosity 
generally outperforms scalar viscosities in preserving geometric symmetries (e.g., sphericity), and suppressing 
viscous heating in homologous spherical contraction. 

The last option we consider is potentially the least diffusive, but also the most unstable or unpredictable 
method due to the matrix inversion operations needed to align the viscosity calculation to the natural 
principal axes of the shock frame. The shock orientation and velocity differences arc determined by the 
eigenvectors R k (S) and eigenvalues X*(S) of the symmetrized strain rate tensor S 1 * = Sym(9jV* — diV % <5*/3). 
In this procedure the artificial viscosity is written as 

Q) = Rl(S) Qf R^Sf , (37) 

where R l j(S) T is the transpose of R l k (S), Qf is the usual von Neumann- Richtmyer viscosity in the shock 
aligned frame where the viscosity tensor has only diagonal elements 

Q\ = I N AV X i (S) (k q Af A* (5) - h C s ) . (38) 
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The shock-aligned length scales Al are computed in a similar fashion, also accounting for local proper frame 
metric distortions, My = g i jAx l Ax : > (assuming no index summation here). Representing the eigenvalues 
and eigenvectors of the matrix My as X l (M) and Rij(M), we define the shock-aligned length scales as the 
diagonal elements of the transformed length scale tensor Al 1 — diag[AZy], where 

Alij = R ia {S) T R ab (M) X bc (M) R cd (M) T R dj (S) , (39) 

and X bc (M) is the diagonal tensor of eigenvalues \ % (M). This expression effectively transforms the principal 
length scales X{M) to the grid axes, then to the principal frame of the strain rate tensor. 

Because the basic strain rate tensor (and especially velocity divergence) does not always distinguish 
between discontinuous shocks and smooth compressive flows, it is sometimes necessary to be more selective 
about where artificial viscosity is applied. As an option we can modify the strain rate tensor as 

diV j => diV j - k L d[-V j , (40) 

where fee is a constant less than unity (typically 0.5), and is a limited gradient. This effectively reduces 
the viscosity levels over smooth flows where the limiter has no affect, while keeping it strongly active over 
shocks where the limited gradient vanishes. In this case we use a hard limiter and set <9f — if any of 
the velocity or velocity gradient components have opposite signs across any of the nearest neighbor zones. 
If both the velocity and its gradient are monotonic in the local neighbor patch, we set d^V^ — diVK 



3.2. NOCD Method 

We take a slightly different and significantly simpler approach in this paper to solving the fully con- 
servative system of equations described in §2.2 than the Ricmann-frec method we implemented in Cosmos 
(Anninos & Fragile 2003). This approach, developed by Kurganov & Tadmor (2000), is also part of the 
NOCD family of solutions, but we find it to be generally less diffusive and less sensitive to Courant restric- 
tions. 

The curvature sources in equations (15) and (16) are updated using finite volume discretization, and a 
first order in time Euler advance method. The divergence terms in (15), (16) and (21) are updated in time 
from u n to u n+1 using a general Ath order method consisting of a sequence of first order Euler steps (Shu 
& Osher 1988). Writing the time update in terms of forward Euler templates, high order time solutions can 
be expressed as 



(i) - 



u 



u 

(m+l) 



= u n + At n S(u n ) 



Vm u n + (1 - Vm )(u^ + At n S(uW)) , 



(41) 



u n+1 = u k 



k — 1 up to k steps, an arbitrary source term S(u), and constant 



(42) 



for ordered sequences m = 1, 2, 
parameters 

J (1/2, -) , 2nd order 
V \(3/4, 1/3) , 3rd order 

At first order, r\ m is not used since that reduces to a simple, single-step forward Euler solver. The time 
update is thus conveniently solved as a sequence of first order Euler cycles, and is easily extended to higher 
order using these simple prescriptions. 
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In semi-discrete form, the single-step solution to a general nonlinear convection equation for a general 
field oj, representing the energy (15), momentum (16), or induction (21) equations, is written in finite volume 
form as 




where F l are the numerical fluxes, a/ = max(aj^, aj) is the maximum local propagation speed, and u>^ = 
lo± + (piuj)± (r 1 ^ — r±) are the limited gradient projections of the cell-centered field io± to the cell faces. 
The subscripts ± refer to either the donor cell center (— ) or the opposite zone (+) across face /, and the 
superscripts ± refer to the face centered projections of the field originating from the corresponding zone. In 
this approach, the transport terms are grouped together with the other fluxes, and the discrete solution (43) 
is applied to all divergence source terms (transport included) simultaneously. 



3.3. Extended AV Method 

The basic internal energy formulation with artificial viscosity can be expanded easily to include an 
additional equation for the conserved energy in the form (15). The general idea in this dual energy approach 
is to extract the thermal component from the total energy, and depending on the accuracy of the result, use 
it to over-write the solution computed directly from the internal energy evolution equation. However, care 
must be used when extracting the internal energy, since it is often the case that adiabats in total energy 
conserving methods can be grossly miscalculated. 

In this scheme, the total energy is used as an option to compute the four-momentum normalization at 
various stages of the solve sequence, as well as the inertia for the artificial viscosity when certain conditions 
are met. Defining £d as the non-thermal or "dynamical" component of the conservative energy, 

DW 2P B W 2 , / 0( , B°B°\ 

£d = ^+ r- +V=9[ 9 00 Pb-— , (44) 



47T 



we write 

(£ - £p)V=g~ W 

rw 2 + (r - i)g 00 (V=g) 

for the internal energy extracted from the conserved energy field, and 



^J^ 2 , (45) 



W(ph + 2P B ) =D + E + W(P +\\Q\\ + 2P B ) 



£ -^99 °(P + Pb) + — V^gB a B°\ (46) 



W V 4tt 

for the inertia and momentum normalization. 

Also, an additional term is added to the total energy evolution equation that accounts for collisional 
dissipation heating arising from the artificial viscosity. The conservative energy equation that we solve takes 
the general form 

d t £ + di (£V l ) + d t (F l ) = S° , (47) 

where the flux F l is now defined as 

1 



F l = V~g ( (g 0j -g°°V*) {(P + P B )5) + Q]) - —{B l B" - B" B"V l ) ) . (4<S) 
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Although we use essentially the same equation for energy (apart from Qj) as the NOCD method described 
in the previous section, we solve it using a different, more conventional operator split approach. Here we 
solve sequentially for curvature effects (d t £ = S°), followed by the transport source (d t £ + di{£V l ) = 0) 
which is updated synchronously with the momentum, internal energy, and magnetic induction advection, 
then finally we apply second order finite volume discretization to evaluate the remaining flux divergence 
term (d t E + d t F l = 0). 

The final step in this extended scheme is to determine whether the total energy solution is known well 

enough to extract accurately the internal energy. We use for a measure of accuracy, the minimal difference 

ratio _ 

'S£\ _ min[E h(T), max(/9e floor W7i(r), E h(T))] N 



£/mi„ max(£) 



N 



(49) 



where 

_ n^+(r-W=g; 

pefloor is the minimum allowable energy density threshold, and the subscript N refers to extending local 
minimum and maximum calculations to include all neighbor zones, adding an extra measure of safety. A 
common problem with total energy schemes in general is that numerical truncation errors can accumulate to 
the point that the sum of different physical contributions often exceeds the total energy (£jj > £), especially 
in kinematic or magnetic field dominated flows, and in the vicinity of strong shocks. This problem is avoided 
by forcing a minimum threshold on either E or E to guarantee positivity, and by preserving the solution 
from the internal energy equation whenever £n > £, or (S£/£ ) m i n < <5 C , where S c is a user defined parameter 
large enough to prevent numerical noise from corrupting the solution. The internal energy is otherwise set 
to E = max(efl oor M /r , E) when (S£/£) min > S c and the extracted thermal component can be trusted. 



4. Code Tests 
4.1. Hydrodynamics 

4.1.1. Shock Tube 

We begin testing with one of the standard problems in fluid dynamics, the shock tube, in which two 
different fluid states are initially separated by a membrane. At t — the membrane is removed and the 
fluid evolves in such a way that five distinct regions appear in the flow: an undisturbed region at each end, 
separated by a rarefaction wave, a contact discontinuity, and a shock wave. Although this problem only 
checks the hydrodynamic elements of the code, as it assumes a flat background metric and ignores magnetic 
fields (magnetosonic shocks will be considered in §4.2), it is still useful for evaluating the shock-capturing 
properties of the different methods. We consider the high boost (W — 3.59, T — 5/3) case from Anninos & 
Fragile (2003). The initial state is specified as pl = 1, Pl = 10 3 , Vj, = to the left of the membrane and 
pn = 1, Pji = 10~ 2 , Vr = to the right. The membrane is located at x = 0.5 on a grid of unit length. 
The results presented here are run using the scalar artificial viscosity with a quadratic viscosity coefficient 
k q = 2.0, linear viscosity coefficient fc; = 0.3, Courant factor fc c /; = 0.3, viscosity multipier — 0, and 
are carried out on fixed, uniform grids of differing resolutions in order to establish the convergence of each 
method. Figures 4, 5, and 6 show spatial profiles of the results at time t = 0.36 on a grid of 800 zones using 
the AV, eAV, and NOCD methods, respectively. Table 1 summarizes the errors in the primitive variables 
p, P, and V for four different grid resolutions (400, 800, 1600, and 3200 zones) and the three different 
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CFD methods (AV, eAV, and NOCD) using the L-l norm (i.e., ||£(a)||i = J2i,j,k Ax i A Vj Az kKj,k - A i,j,k\> 
where a" - k and A^. k are the numerical and exact solutions, respectively, and for one-dimensional problems 
the orthogonal grid spacings are set to unity). All three methods converge at approximately first order as 
expected for this class of problem, and the errors are consistent with those reported for our previous code 
(Anninos & Fragile 2003). 



4-1-2. Shock Collision 

A second test presented here is the wall shock problem involving the shock heating of cold fluid hitting 
a wall at the left boundary (x — 0) of a unit grid domain. The initial data are set up to be uniform across 
the grid with adiabatic index V = 4/3, pre-shocked density p\ = 1, prc-shocked pressure Pi = (T— 1) x 10~ 8 , 
and velocity V\ — ~v init . When the fluid hits the wall a shock forms and travels to the right, separating the 
initial pre-shocked conditions from the post-shocked state (p2, P2, V2) with solution in the wall frame 

VS = T7T-, 51 

p 2 - p\Wx 

i , 2 = p 2 (r-i)(Wi-i), (52) 



P2 = Pi 



r + 1 ' r o*-i) 



(53) 



r- 1 r- 1 

where Vs is the velocity of the shock front, and the pre-shocked energy and post-shocked velocity are both 
assumed negligible (ei, V2) — > 0. All of the results in this section are performed on a 200 zone uniformly- 
spaced mesh and run to a final time of t = 2.0. For the AV and eAV methods, we use the scalar viscosity 
with k q = 2.0 and fc; = 0.7. The Courant factor is set to k c fi — 0.3 for all methods. 

Figure 7 plots the mean-relative errors in density, which are generally greater than errors in either the 
pressure or velocity, as a function of boost factor. Although we are not able to extend the AV method 
reliably (which we define by a 10% mean error threshold and increased sensitivity to viscosity parameters) 
beyond Uj n , t ~ 0.95, the eAV and NOCD methods are substantially more robust. As shown in Figure 8, 
both methods can be run up to arbitrarily high boost factors (vi n n > 0.99999) with mean relative errors 
typically remaining below 2% with no significant increasing trend. As noted previously (Anninos & Fragile 
2003), the errors for the AV method can be improved significantly by either lowering the Courant factor 
or increasing the viscosity coefficients. However, the sensitive dependence on these parameters detracts 
considerably from the attractiveness of this method for this class of problem. The errors for the eAV method 
can also be improved by adjusting the viscosity coefficients, although we find that it performs well even with 
the "standard" values used here. 



4-1.3. Boosted Shock Collision 

An elaboration of the wall shock problem from the previous section is the collision of two boosted 
fluids. The fluids flow toward each other and collide, each shocking against the other and forming a contact 
discontinuity. In the center-of-momentum frame, in which the contact discontinuity between the fluids 
is stationary, this system is equivalent to a pair of opposing wall shocks, each impinging on the other. 
However, by boosting this system's center-of-momentum frame with respect to the grid frame, one devises a 
very rigorous test of the Lorentz invariance of the code under nonsymmetric conditions, with multiple jump 
discontinuities and highly relativistic shock velocities. 
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In the center-of-momentum (primed) frame, a fluid with initial state (p±, Pi, V{ > 0) flows from the 
left boundary while a fluid of state (p 4 , P4, V 4 ' < 0) flows from the right; we assume each fluid is cold: 
Pi = P 4 = 0. Upon collision, fluid 1 is shocked into a state (p2, P2, V2) while fluid 4 is shocked into state 
(p3, P 3 , U 3 '), where the fluids are numbered sequentially from left to right. In this frame, each shocked fluid 
comes to rest; V^ = V£ = 0. This implies force equilibrium between these fluids; P2 = P3. Using equations 
(52) and (53), pressure balance gives 

Pl (TW{ + 1){W{ - 1) = Pi (TWi + l)(Wi - 1) . (54) 

The proper densities of the shocked fuids (from eqn. 53) are 

rwj + 1 .... 

P2 = Pi—^ — t- (55) 



r- 1 

P3 = Pa r _ 1 , (56) 



and the specific energies are (e.g. Hawley et al. 1984) 

e 2 = W[-l (57) 
e 3 = Wi - 1 . (58) 

Pressure balance again implies p2&2 = P3£3- Also, from equation (51) the reverse shock velocity, U r ', between 
materials 1 and 2, and the forward shock velocity, Vj, between fluids 3 and 4 are expressed as 

(iu{-i)(r-i) 

y " = (59) 

Vi = - ( ^- 1)(r - 1} (60) 
' 1U^U 4 ' 1 ; 

Typically the grid (lab frame) velocities V\ and V4 are known and the velocity of the center-of-momentum 
frame and contact discontinuity, Vd, must be solved for numerically by substituting the boost transformations 

W[ = W d Wi(l - ViV d ) , Wi - W d Wi(l + V 4 V d ) (61) 

into equation (54), where Wi = (1 — V- 2 )^ 1 ^ 2 . The system can also be boosted into a desired frame by 
standard velocity addition. 

Figures 9, 10 & 11 compare the analytic and numerical solutions, using the eAV method, for three 
examples of boosted collisions. We find very good agreement to the solutions (Table 2) when the opposing 
primed- frame momenta of each fluid is symmetric (Figure 9), asymmetric (Figure 10) and for a more extreme 
relativistic case with Lorentz factors of about 100 (Figure 11). Fractional errors for density in the shocked 
region are typically Ap/p <~ 10~ 4 while both proper energy and boost factor fractional errors are of order 
10~ 5 . The tail in mass density on the downwind (right hand) side of the shocked region is a small proportion 
(3% for Figures 9 & 10 and 1% for Figure 11) of the total shocked mass and converges to zero with increased 
zoning. The NOCD method also gives good agreement, but tends to be slightly more diffusive, with a larger 
mass tail on the trailing shock. In the most extreme relativistic case, shown in Figure 11, the numerical 
results are slightly behind (to the right of) the analytical solution, corresponding to temporal delay in the 
onset of the numerical shock of At w 10~ 4 . These highly relativistic shocks require very fine zoning and 
thus AMR is extensively used in order that these runs be computationally feasible. Up to 12 levels of 
refinement are used in the highest boost test. In these examples a simple density threshold is used, above 
which refinement is triggered; however, more complex (and efficient) triggering methods such as discussed 
in the next section have also been successfully employed. 
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4-1. 4- Blast Waves 

A further application of the ultra-relativistic capabilities of the eAV and NOCD methods demonstrated 
in the previous section and one that is of particular astrophysical significance is the relativistic blastwave. 
Much like its Newtonian counterpart, the Scdov- Taylor blastwave (e.g. Landau & Lifshitz 1959), there is 
a self-similar analytical solution for the ultra-relativistic blastwave, first described by Blandford & McKee 
(1976). 

The analytical solution of Blandford & McKee (1976) depends on the initial total energy in the blastwave, 
Ebm, the Lorentz factor with which its shock front initially expands, 7bm, and the ambient density, pbm, 
into which it expands. From these one can define an initial radius of the blastwave 

( UE BM \ 1/3 

TBM = o 2 • 62 

The solution for radii, < r < r BM , is based on the similarity variable 

X (r) = 1 + 8 7bm (1 - r/r BM ) ■ (63) 

The coordinate density and energy, as used by Cosmos++ , and the Lorentz factor are 

D{r) = 2V^Pbm7ImX(0" 7/4 (64) 
E(r) = V ^ P bm7Im (V27bmx(0^ 23/12 - 2 X (Q- 7/4 ) (65) 

7(0 = f^jfa (66) 

where we have interpreted equation (29) of Blandford & McKee (1976) as the radial component of the 
4-velocity in order that 7 > 1. For all r > tbm, D(r) = ^/~ r gpBM, 7 (v) = 1 and the energy is set to a 
numerically insignificant value, typically E(r) = 10~ 4 D(r). 

The relativistic blastwave is characterized by a very thin, Ar oc t"bm/(87 bm ) [eqn. 63], shell of matter 
and energy rapidly expanding into an external medium. Very fine zoning is required to resolve the shell, while 
relatively few zones are necessary in the vacuous bubble enclosed by the shell or the external medium. Thus 
this problem is ideally suited for adaptive mesh refinement. For example, a blastwave with an initial Lorentz 
factor of 7bm = 10 would require 800 zones across its radius just to ensure that the shell is represented by 
a single zone. Because of the steep gradients behind the shock, to resolve the shell requires / ros ~ 100 times 
as many zones, which becomes a large computational problem. This is compounded by the fact that typical 
simulations evolve the blastwave over distances cx tbm, thus requiring about r^yi/Ar <~ t-qmI^cs/ [k c fi At) ~ 
800f rcs /k c fi ~ 10 6 time steps, where fc c // is the Courant factor. For 7bm > 10, stability typically requires 
k c f 1 ~ 0.1. Therefore adaptive mesh refinement, with zones concentrated into the shell, is a practical 
necessity for relativistic blastwaves. 

The relativistic blastwave, with its thin shell, steep gradients, and sharp cusp at the shock front, is 
challenging to simulate. The method developed for these problems is to refine a zone if any of three criteria 
are met: First, all zones for which the density p > /thresh/^max are refined to the fullest allowable extent, 
where typically /thresh = 0.9. This condition is required to keep the cusp at the shock front as resolved as 
possible. Second, the normalized derivative of the proper density field, (Ax /7j)V p, must remain below a 
threshold, where Ax is the zone size and p is the average value. This is effective at maintaining the steep 
gradient behind the shock. Finally, the normalized curvature of this field, Ax\ V 2 p|/| Vp|, must remain below 
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a maximum threshold. This condition puts zones ahead of the advancing shock due to the high curvature 
at the shock discontinuity, and also smooths the regions tagged by slope. If p < /thrcshPmax and both 
normalized slope and curvature of a zone fall below their respective thresholds, then that zone is tagged for 
derefinement. 

As shown in Figure 12, Cosmos++ is able to evolve a blastwave with -Ebm = 10 51 ergs from an initial 
Lorentz factor of 7bm — 30 to a non-relativistic blastwave. Such a simulation is demanding, requiring ~ 10 6 
cycles to evolve the relativistic phase. This simulation is evolved on a base mesh of 100 zones with initially 
17 levels of allowed refinement. The number of allowed refinement levels is periodically decremented nine 
times as the blastwave decelerates and broadens. Thus only eight levels of refinement are employed during 
the non-relativistic phase. As a result, the time step during the final non-relativistic phase is 2 10 w 1000 
times larger than the initial timestep. 

Figure 13 shows the self-similarity of the relativistic and non-relativistic phases of the blastwave evolution 
by scaling and super-imposing several density profiles for each phase. One can see excellent agreement with 
the analytical solution in both cases. For the relativistic phase, the eAV method (shown in the left plot) 
tends to more robustly evolve the analytical solution, with less sensitivity to time step or AMR refinement 
prescriptions. However, the NOCD method (shown in the right plot) more reliably transitions to the non- 
relativistic limit; for instance capturing the shock discontinuity in density with ~ 1 % accuracy as opposed to 
<~ 10 % accuracy for the eAV method. This might be due to the enforced energy conservation of the NOCD 
method. The eAV method typically only loses 1 - 2% of the total energy over the course of the evolution, 
but this may occur primarily in the shocked shell, thus causing a slight drift in the solution. 



4.2. Magnetohydrodynamics 

4-2.1. Alfven Wave Propagation 

The class of linear Alfven waves described by De Villicrs & Hawley (2003a) provide an excellent test 
of the method of characteristics for magnetic fields subject to transverse or shearing mode perturbations. 
Considering a general wave function fix — v A t), solutions to the linear perturbation equation with a fixed 
background field B x and constant velocity V x in Minkowski spacetime yields for the transverse components 

yy{x,t) = ( 1 —^.)f( x -v A t) + ( 1 -±^.)f( x -v+t) 1 (67) 



B»{x,t) = ^[f(x-v A t)-f(x-v+t)] , (68) 



with parameters 



and Alfven speed 



_ y*±^r? + W -2 _^ / ||i?|| 2 _ / 2(r-l) £ 

I + 77 2 y47r^+||B|| 2 y + re) + 2(r - l)e " 1 ' 

We consider two cases of linear Alfven waves: a stationary background (V x — 0, case A) in which the pulses 
travel in opposite directions with equal amplitudes, and a moving background (V x — 0.1c, case B) where 
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the pulses split into asymmetrical waves. These two cases correspond to Models ALF1 and ALF3 of De 
Villicrs & Hawley (2003a). In both cases, the fluid is initialized with a uniform density p = 1, specific energy 
e = 10~ 2 , and adiabatic index V — 5/3. The transverse magnetic field components (B v = B z ) are initially 
zero. The longitudinal magnetic field component (B x ) is set by our choice of (3 (/3 = 0.001 for case A and 
0.01 for case B). The transverse velocity is initialized with a square pulse: V v = 10~ 3 for 1 < x < 1.5, 
V y = — 10~ 3 for 1.5 < x < 2, and zero elsewhere. The problem is run on a grid 3 units in length. 

For case A, the Alfven velocities are |t>^| = 0.963 (W = 3.76); for case B, v\ = 0.792 (W+ = 1.64) and 
= —0.705 (W~ = 1.41). Figure 14 shows the numerical results for the artificial viscosity (AV) method on 
a grid with 1024 zones overlaid with the analytic solution. The eAV method gives identical results for this 
test and is not shown. The errors in the plateaus for the stationary background case (case A) are extremely 
small (< 0.005%); for the moving background (case B), the errors are somewhat larger, though still quite 
small (< 0.1%) 

4-2.2. Magnetosonic Shock Tube 

Next we perform another set of numerical shock tube tests, this time including magnetic fields. Similar 
to the hydrodynamic shock tube test, these problems combine strong shocks and rarefaction features, thus 
fully stressing the code in the limit of flat space. We initialize two versions of this test, both from Komissarov 
(1999). In the first test, the initial state is specified as p L = 1, P L = 10 3 , V L = 0, B X L = 1, B\ = B Z L = to 
the left of the membrane and p R = 0.1, Pr = 1, Vr = 0, B R = 1, B V R = B R = to the right. Note that since 
the magnetic field in this case is parallel to the flow and continuous across the discontinuity, it should not 
play a dynamical role. In this sense, the problem is really a hydrodynamical shock tube similar to the one 
considered in §4.1.1. Nevertheless, this test appears frequently in the literature and is worth reproducing here 
to confirm consistency in the MHD solver. For the second test, the initial state is Pl = 1, Pl = 30, Vl = 0, 
B\ = 20, B X L = B Z L = to the left of the membrane and p R = 0.1, P R = 1, V R = 0, B X R = B y R = B Z R = 
to the right. Here, the initial discontinuity of the field allows it to play a dynamical role in the rarefaction 
and contact wave regions. These tests are run using the scalar artificial viscosity with a quadratic viscosity 
coefficient k q — 2.0, Courant factor k c fi — 0.3, and are evolved on fixed, uniform grids. Case 1 uses a linear 
viscosity coefficient ki — 0.3, whereas case 2 does not use the linear term (fc; = 0). We do not calculate the 
analytic solutions for these tests, although a direct comparison can be made between our Figures 15, 16, and 
17, and Figure 6 of Komissarov (1999) since we use the same resolution (400 and 500 zones for shock tubes 
1 and 2, respectively) as that work. We point out that our different methodologies all give similar results 
as the Godunov-type scheme used in Komissarov (1999), showing many of the same pathologies (the small 
post-shock features in shock tube 1 and the kink at the rarefaction edge in shock tube 2), although our AV 
and eAV methods appear to do better at capturing the density plateau in the Lorentz-contracted shell of 
material immediately behind the shock. 

4.2.3. Bondi Flow 

As a test of hydrodynamic and MHD flows in spacetimes with nontrivial curvature, we first consider 
radial accretion onto a compact, strongly gravitating object, in this case a Schwarzschild black hole. The 
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analytic solution (Michel 1972) is characterized by a critical point in the flow as 

K) 2 = M/2r c , (71) 

v; - KftH-aiO^ ,;^,;^ : (72) 

where u r is the radial component of the fluid 4- velocity, V c is the sound speed at the critical point, M is the 
mass of the black hole, n = 1/iY — 1) is the polytropic index, and T = P/p. The remainder of the flow is 
described through the equations 



T n u r r 2 = Ci , (73) 



[l + (l + n)T] 



r 



= C 2 . (74) 



Following Hawlcy et al. (1984), we fix C\ and C2 by choosing the critical radius of the solution r c — 8GM/c 2 
and fixing T = 4/3. We also fix the density at the critical radius (p c ) such that M = \-Kr 2 p c u r c = —1. 

The physical domain of our simulation extends from r = 0.98rsH to r = 20GM/c 2 , where tbh = 
2GM/c 2 is the radius of the black-hole horizon. Here we use Kerr-Schild coordinates, which allow us to place 
the inner boundary of the computational domain inside the horizon. The radial coordinate is replaced by a 
logarithmic coordinate 77 = 1 + ln(r/r^ff), and the problem is evolved over a time interval At = 100GM/c 3 . 
We measure the convergence of our solution at four different resolutions (32, 64, 128, and 256 zones) using 
the one-dimensional L-l norm of p. For 256 zones, the respective errors for the AV, eAV, and NOCD methods 
are 1.14 x 10~ 3 , 1.08 x 10~ 3 , and 4.66 x 10~ 4 . The convergence rates are between first and second order: 
1.3, 1.4, and 1.4 for the AV, eAV, and NOCD methods, respectively. These rates reflect a slight degradation 
of truncation order as expected for curvilinear meshes since we do not currently construct volume-centered 
centroids nor high order face-centered interpolants that affect gradient operators, for example. Nevertheless, 
these convergence rates are generally consistent with the convergence reported by De Villiers & Hawley 
(2003a). 

We can extend this test to the GRMHD regime by adding a radial magnetic field. Inclusion of such 
a field (satisfying d r B r = 0) does not alter the analytic solution for any of the primitive fields (p, P, or 
u r ). We should point out, however, that this treatment does not satisfy the whole set of Maxwell equations 
(Anton et al. 2005), yet it serves as a non-trivial numerical test of the magnetic field terms in the code and 
is well documented in the literature. The magnitude of the magnetic field is set by ||_B 2 ||/(47rp) = 10.56 
(/? s=s 1) at r = r c . With magnetic fields included, the magnitudes of the errors increase to 3.89 x 10~ 3 and 
1.36 x 10~ 3 for the AV and eAV methods with 256 zones, but the convergence rates remain the same - 1.3 
and 1.4, respectively. We note that the divergence error in the magnetic field does not build up appreciably 
during the course of these runs - \diB l \ < 10~ 14 in all cases. Since we intend primarily to use either the 
AV or eAV methods for this type of research application, we have not extended the NOCD algorithm to 
account for the proper characteristic magnetohydrodynamic speed in arbitrarily curved spacetimes, so wc 
do not report results from that method in this or the following section. 



4-2.4- Magnetized Black Hole Torus 

Next we consider the astrophysically interesting problem of a magnetized torus of gas orbiting a rotating 
black hole. Although there is no known analytic solution for this problem, it is sufficiently well documented 
in the literature (e.g. De Villiers & Hawley 2003a; Gammie et al. 2003; Anton et al. 2005) to serve as a useful 
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test of the code. It also represents one class of problems for which Cosmos++ is intended to be used. We set 
this problem up using the same parameters as model SFP-2D of Dc Villicrs & Hawley (2003b) ; specifically, 
the spin of the black hole is a/M — 0.9, the specific angular momentum of the torus is l/M = 4.3, the surface 
potential of the torus is (u t )i n = —0.98, and the equation of state of the gas is set by the polytropic constant 
k = 0.01 and index T = 5/3. The initial setup of the torus is illustrated in Figure 18(a), which shows a plot 
of the logarithm of p. 

To this torus we add initially weak poloidal magnetic field loops to seed the magneto-rotational instability 
(MRI, Balbus & Hawley 1991). The initial magnetic field vector potential is (De Villiers & Hawley 2003a) 



The non-zero spatial magnetic field components are then B r = —deA^ and B = drA^. These poloidal field 
loops coincide with the isodensity contours of the torus. The parameter p cut — 0.5p max is used to keep 
the field a suitable distance inside the surface of the torus. Using the constant k in equation (75), the field 
is normalized such that initially — P/(\ \B\ | 2 /87r) > 0q — 2 throughout the torus. This initialization is 
slightly different than De Villiers & Hawley (2003b), who use a volume integrated to set the field strength; 
the difference is such that 0o = 100 in their work is roughly comparable to 0o = 2 here. 

In the background region not specified by the torus solution, we initially set up a cold (e = 10~ 6 e max ), 
low-density (p = 10 _6 p max ), static (V 1 = 0) gas, where e max and p max are the initial internal energy and 
mass densities at the pressure maximum of the torus (at r center — l5.3GM/c 2 ). These values for e and p 
are then used as floors throughout the simulation. Obviously once the evolution begins, the background 
is no longer static. At first it falls toward the hole, but ultimately most of the background region is filled 
with a magnetized "corona" and a low density, high outflow launched by pressure forces and the opening 
of magnetic field lines. Because the background gas always remains very low density, it does not have a 
significant dynamical effect on the torus, although its role in real astrophysical systems is still not well 
understood. 

For this test, we perform a two-dimensional, axisymmetric simulation using our AV method. We restrict 
this test to a single level fixed mesh with a resolution of 128 2 . In later applications, we plan to employ 
the full three-dimensional, refined grid capabilities of Cosmos++ . This simulation uses a logarithmic 
radial coordinate of the form T] = 1 + ln(r/rs//) and a concentrated latitude coordinate X2 of the form 
= X2 + |(1 — h) sin(2x2) with h = 0.5. The grid covers the angular scale 0.027T < 9 < 0.987T and has radial 
boundaries at r m i n ~ Q.Q&tbh and r max — 120M. We use outflow boundaries at both locations (radial fluid 
motion onto the grid is not allowed). The chosen grid resolution gives a zone spacing of Ar s=a 0.05GM/c 2 
near the inner radial boundary and Ar m 0.5GM/c 2 near the initial pressure maximum of the torus. This 
can be compared with the characteristic wavelength of the MRI, Xmri = 2ttva/Q ~ 2.5GM/c 2 near the 
initial pressure maximum. 

Figure 18 shows the mass density distribution at t/(GM/c 3 ) = 0, 1250, and 2300 (roughly 0, 3, and 
6 orbits at the initial pressure maximum). For a more quantitative comparison, Figure 19 shows the time 
history of the mass accretion rate through the inner radial boundary of the grid. Comparing with Figure 18 
of De Villiers & Hawley (2003b), we see similar amplitudes and frequencies of variability in the accretion flow, 
particularly for the first 4 orbits. We also find a similar value for the total mass accretion: AM/Mq = 0.13 
after 6 orbits here versus 0.14 after 10 orbits in Dc Villicrs & Hawley (2003b). The accretion tapers off to 
a low value after about 3.5 orbits as expected since MRI turbulence cannot be sustained indefinitely in an 
axisymmetric simulation, due to Cowling's antidynamo theorem. Throughout this simulation the divergence 
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of the magnetic field is held below a tolerable limit IdiB 1 ] ~ 0.01) through the use of the divergence 
cleanser. Importantly, this divergence error does not increase with time. 

5. Discussion 

We have developed a new object oriented code Cosmos++ for solving the fully general relativistic magne- 
tohydrodynamics equations on adaptive, unstructured meshes using discrete finite volume and dimensionally 
unsplit methods. Three basic numerical schemes have been implemented and tested using both internal and 
total conservative energy formulations: The first evolves internal energy with artificial viscosity for shock 
capturing; the second uses a nonoscillatory central difference scheme to solve the fully conservative (energy 
and momentum) form of equations; and the third approach combines the internal and conservative energy 
equations with artificial viscosity methods to achieve greater accuracy in highly relativistic regimes. 

We find by comparing the different methods presented here with other published results, including 
Riemann-solver based codes, that, despite their simplicity, artificial viscosity methods perform quite well for 
low to moderately boosted flows, (V/c ~ 0.95 in the shock tube and wall shock tests). This is consistent 
with our conclusions from earlier work using dimensionally split, finite difference methods on structured grids 
(Anninos & Fragile 2003). However, it is well known that for higher velocity flows, traditional AV methods 
tend to break down. We have demonstrated that the basic AV approach can be extended easily to the ultra- 
relativistic regime by simply incorporating a dual energy formalism (the eAV method) to solve both internal 
and conservative energy equations with standard operator split procedures. This eAV procedure results 
in significantly improved shock and wave capturing capabilities that allows AV schemes to model flows at 
arbitrarily high boosts. We have presented stable, accurate solutions for strong shock collision interactions 
with boosts (velocities) easily exceeding 200 (0.99999c). 

The NOCD method implemented here, although of a slightly different family of algorithms than wc 
adopted in our previous code, also works quite well in the high boost regime. In fact this method has signif- 
icantly less numerical diffusion and is less sensitive to Courant restrictions than our earlier implementation 
(Anninos & Fragile 2003). The eAV and NOCD methods thus provide robust alternatives to simulating 
highly relativistic flows since they are comparable in accuracy, over the entire range of velocities we have 
simulated, to more complicated Ricmann-based codes. A similar conclusion regarding central difference 
schemes was reached by Anninos & Fragile (2003) and Lucas-Serrano et al. (2004). Here we have extended 
the scope of tests and validity of this class of methods to relativistic magnetohydrodynamics. 
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grant AST 0307657. Support was also provided by the National Science Foundation under the following NSF 
programs: Partnerships for Advanced Computational Infrastructure, Distributed Terascale Facility (DTF) 
and Terascale Extensions: Enhancements to the Extensible Terascale Facility. 
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*N 2 

Fig. 2. — Compilation of the local attributes stored in each cell object to specify the geometry and connectiv- 
ity data. Each cell carries dynamically allocated pointers to its parent cell (*p), four (eight) children (*C n ), 
and four (six) neighbors (*N n ) in 2D (3D). In order to accommodate arbitrary cell shapes and geometries, 
the node positions x l n , face area normals A l n (centered at each cell face and directed toward the cell center), 
and cell volumes are also stored and allowed to change during the evolution. The child and neighbor lists 
are ordered as shown in all the zone objects for fast and easy memory access. This data structure allows for 
non-orthogonal, arbitrarily distorted and unstructured evolving meshes. 
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Fig. 3. — Schematic of the sub-zonal decomposition used in constructing upwind (U) and downwind (D) 
stencil geometries for enforcing monotonicity in the gradient extrapolants. In this example, the two shaded 
triangular domains comprise the total upwind polygonal sub-zone as determined by the sign of the inner 
product of the face area normal and the face-centered velocity vectors {AiV 1 > 0). The two unfilled triangular 
domains separated by a dotted line represent the composite downwind polygonal. Gradients are computed 
separately on the three different cell blocks (upwind, downwind, and average which includes the entire 
original quadrilateral cell) using finite volume contour integrals with fields appropriately averaged to each 
of the cell faces. The final gradient is a combination of the three sub-zonal constructs and limited according 
to the procedures discussed in the text for monotonicity. In 3D, the sub-zone domains are constructed from 
tctrahcdral elements. 
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Fig. 4. — Results at time t = 0.36 for the shock tube test using artificial viscosity (AV) and 800 zones. The 
data points in this plot have been sampled to reduce overcrowding. Only 400 points are shown. The solid 
line shows the analytic solution. 



- 25 - 




- 26 - 




Fig. 6. — As Figure 4 but with the non-oscillatory central difference (NOCD) scheme. 
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Fig. 7. — Mean relative errors in density for the AV, eAV, and NOCD methods as a function of boost for 
the relativistic wall shock problem. All calculations were run using 200 zones up to time t = 2.0. The AV 
and eAV results can be improved significantly and brought closer in alignment with the NOCD results by 
reducing the Courant factor or increasing the viscosity coefficients over the canonical values we have chosen. 
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Fig. 8. - Density plots for different infall velocities in the wall shock test using the eAV (squares) and NOCD 
(triangles) methods. Results for the AV method (circles) are only included for V = —0.9. The resolution is 
200 zones, the displayed time is t = 2.0, and the graphs are zoomed in at the shock fronts to distinguish the 
different solutions. 
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Fig. 9. — Two ideal fluids, each with initial proper density p = 1, proper energy pe = 1CP 8 , and T = 5/3, 
collide initially at x = 0.05. The fluids move in opposite directions, each with W = 5 in the center-of- 
momentum frame. The observer is boosted to the right at W — 3 so the fluid moves at speeds up to 
V/c ~ 0.999, and the shocked fluid region can be seen to move to the left over a sequence of five times, t = 
0.0, 0.01, 0.02, 0.03, 0.04. The numerical results are solid lines and show good agreement with the analytical 
solution (Tabic 2) shown in dot-dashed lines. The ticks at the top of the plot are zone positions for the last 
time snapshot. This problem uses a base resolution of 60 zones, with eight allowed levels of mesh refinement. 
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Fig. 10. — As in Figure 9 except the two ideal fluids have unmatched initial proper densities p\ = 2, pi = 1. 
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Fig. 11. — As in Figure 9 except here W = 10 and the observer is boosted to the right at W = 5. The base 
resolution is 60 zones and 12 levels of mesh refinement are allowed. 
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r (lightyears) 

Fig. 12. — A spherical relativistic blastwave with energy Ebm = 1 foe = 10 51 ergs, 7bm = 30, plows into an 
external density of pbm = 1 baryon cm~ 3 = 1.27 foe lightyear^ 3 from an initial radius tbm = 0.084 lightyear. 
The material immediately behind the shock is initially shocked to p = 2 3 / 2 7bmPbm, 7 = 7bm/\/2 and evolves 
as p oc 7 — 1 oc r~ 3 / 2 . Once the blastwave decelerates to 7 ~ 1 it transitions to the non-relativistic Sedov- 
Taylor blastwave for which the shocked material evolves as p = pbm (T + l)/(r — 1) = 8.9 foe lightyear -3 , 
T = 4/3, and 7— 1 oc v 2 oc r~ 3 (Landau & Lifshitz 1959). This simulation uses the NOCD method. For this 
figure the dump frequency was halved periodically to allow late time profiles to be resolved. 
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Fig. 13. — Self-similarity of the relativistic (left) and non-relativistic (right) phases of the blastwave of 
Figure 12. As discussed in the text, the relativistic (left) panel demonstrates the eAV method, while the 
non-relativistic (right) panel uses the NOCD method. At left, the initial relativistic proper density profile, 
p = D/W, (dot-dashed line) of eqn. (66) is plotted with four subsequent profiles at ~ 100, 000 cycle intervals. 
Each profile is scaled in radius to align with the initial shock front and the density is scaled by p oc r~ 3 / 2 . 
Save an initial decrease in the magnitude of the peak of about 10%, and a gradual increase in the density 
of the lower velocity tail, the self-similar profile is well maintained. At right, the Scdov solution (dot-dashed 
line, e.g. Landau & Lifshitz 1959) with T = 4/3 is plotted with six density profiles from the non-relativistic 
phase of the run from Figure 12. The radii, ranging from 1 to 1.2 lightyears for times 1.33 to 2 years 
respectively, are scaled to the radius of one lightyear. The shocked density p — pbm(X + l)/(r — 1) = 8.9 foe 
lightyear -3 is not scaled, but is naturally captured by the code to within a few percent. The scaled profiles 
are virtually identical to each other and well reproduce the Sedov profile. Tick marks at the top of the plots 
show the mesh nodes for a typical profile. The left (right) plot has 17 (8) allowed levels of refinement. 




Fig. 14. — Tranverse velocity V y for two cases of the Alfven wave test: (a) case A at t = 0.9 and (6) case B 
at t = 1.1. The numerical resolution is 1024 zones; the solid lines are the analytic solutions. 
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Fig. 15. — Left panel: Magnetosonic shock tube problem 1 from Komissarov (1999). The magnetic field 
is normal to the initial discontinuity (at x = 0) so it does not play a dynamical role. The symbols show 
the numerical solution at t = 1 for 400 zones using the AV method. Right panel: Magnetosonic shock 
tube problem 2 from Komissarov (1999). The initial magnetic field to the left of x — is parallel to the 
discontinuity while the field is absent to the right. The symbols show the numerical solution at t — 1 for 
500 zones using the AV method. 
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Fig. 17.— 



As Figure 15 but with the non-oscillatory central difference (NOCD) scheme. 
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Fig. 18.— Plot of gas density p at (a) t = 0, (6) t = 1250M = 3r or6 , and (c) t = 2300Af = Qr orb for 
the Kerr-Schild form of the metric. The density is scaled logarithmically over 4 orders of magnitude and 
maintains the same scale in all three panels. 
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Fig. 19. — Mass accretion rate M at the inner radial boundary r m i n normalized by the initial mass of the 
torus. 
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Table 1: L-l Norm errors in density, pressure, and velocity for the hydrodynamic shock-tube test at time 
t = 0.36. Convergence is approximately linear as expected for problems with shock discontinuities. 
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1.0 


1.1436 
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14.0 


3.0 
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14.0 


1.741 
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28.86 


2.0 


1.1436 
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24.28 
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16.20 
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28.86 


1.0 


1.2558 


9.973 


26.5 


5.0 


238.5 


26.5 


2.565 


1.0 


98.74 



Table 2: Analytical results for boosted collisions shown in Figs 9, 10 & 11 (see §4.1.3). Initial conditions are 
given by pi, W\, p4, W4. Wd is solved for by Equations (54) & (61) in the center-of-momentum frame, and 
then boosted into the lab frame. Shocked proper densities pi and pz are given by eqns. (56), and the proper 
energy density of the contact discontinuity is = £2^2 = £3^3 with eqns. (58). The boosts of the reverse 
and forward shock fronts, W r and Wf, are calculated from their center-of-momentum velocities (eqn. (60)). 



